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I  SUMMARY 

The  research  performed  under  the  contract  during  the  period  1  May 
through  31  October  1981  can  be  divided  Into  three  main  topics,  coupling 
of  surface  waves  in  laterally  inhomogeneous  source  regions  to 
teleseismic  propagation  paths,  using  regional  waveforms  to  determine  the 
source  parameters  of  moderate-size  earthquakes  and  the  application  of 
the  Kir choff -Helmholtz  integral  to  seismic  problems. 

In  Section  II,  the  Representation  Theorem  (RT)  is  used  to 
numerically  evaluate  the  effectiveness  of  two  commonly  used  algorithms 
for  modeling  Rayleigh  wave  propagation  across  lateral  inhomogenities. 
Of  the  two,  the  conservation  of  lateral  energy  flux  approximation  most 
closely  matches  the  maximum  peak-to-peak  amplitudes  as  seen  through  an 
LP-LRSM  instrument.  Another  interesting  result  is  that  in  using  the  RT 
to  model  seismic  events,  the  short-periods  are  dominated  by  the 
displacement  forcing  functions  and  the  long  periods  are  controlled  by 
the  stress  forcing  functions.  This  is  in  agreement  with  other 
investigators  in  that  they  found  that  body  waves  were  controlled  by  the 
displacement  forcing  functions  and  surface  waves  by  the  stress  forcing 
functions. 

In  Section  III,  a  procedure  for  the  systematic  determination  of 
source  parameters  from  regional  body  waves  is  presented.  A  least 
squares  inversion  based  on  a  cross-correlation  of  the  data  and 
synthetics  Is  used  to  determine  the  fault  mechanisms  of  a  profile  of  the 
Pn£  synthetics  and  five  earthquakes.  The  synthetics  are  for  a  western 
U.S.  model  which  seems  to  be  more  than  adequate  for  most  continental 
regions.  Three  of  the  earthquakes  are  in  the  western  U.S.  Two  other 
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earthquakes,  one  in  Baffin  Bay  in  the  Arctic  and  the  other  in  Turkey 
have  been  included  to  demonstrate  that  the  Green's  functions  are  not 
unique  to  the  U.S.  Both  dip-slip  and  strike-slip  events  are  determined. 
It  is  shown  that  the  inversion  parameters  are  fairly  insensitive  to 
small  changes  in  crustal  thickness,  Pn  velocity  and  mean  crustal 
velocity.  The  inversion  precedure  only  requires  a  small  data  set  and  is 
particularly  ideal  for  strike-slip  earthquakes. 

A  numerical  method  for  evaluating  the  Kirchof f-Helmholtz  Integral 
with  application  to  seismic  problems  is  presented  in  Section  IV.  The 
method  is  applied  to  the  calculation  of  reflections  from  mountains 
topography  and  the  modeling  of  spall  associated  with  nuclear  events. 
For  the  latter,  the  calculations  produce  travel-time  and  amplitude 
anomalies  consistent  with  observations  from  underground  nuclear  blasts. 
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II 

EVALUATION  OF  MIXED  PATH  RAYLEIGH  WAVE  TECHNIQUES  USING  THE 
REPRESENATION  THEOREM  I.  VERTICAL  POINT-FORCE  SOURCES. 


by 

Peter  Glover  and  David  G.  Harkrider 


INTRODUCTION. 

In  our  semi-annual  report  for  May  1981  we  presented  some 
preliminary  results  comparing  a  Representation  Theorem  approximation 
with  two  previously  used  analytical  techniques  for  calculating  Rayleigh 
waves  generated  by  explosions  in  mixed-path  media. .  The  Representation 
Theorem  method  that  we  presented  was  approximate  in  that  we  used  a  whole 
space  formulation  for  the  forcing  functions  and  coupled  them  to 
halfspace  Green's  functions.  Thus  the  interactions  of  the  P-waves 
generated  by  the  explosion  with  the  free  surface  within  the  source 
region  and  the  contact  between  the  paths  were  ignored.  The 
discrepancies  that  we  found  between  the  Approximate  Representation 
Theorem  (ART)  results  and  the  two  analytic  methods,  referred  to  as  the 
Conservation  of  Lateral  Energy  Flux  (CLEF)  and  Unit  Transmission 
Coefficient  (UTC)  methods  respectively,  together  with  the  large 
differences  in  amplitude  between  the  CLEF  and  UTC  results,  indicated 
that  we  need  to  make  the  comparison  of  the  two  methods  using  the  exact 
Representation  Theorem  (RT)  results.  To  do  this  we  require  forcing 
functions  for  an  explosion  in  a  half-space  which  we  propose  to  model 
using  the  finite-element  code  SWIS  (Frazier  and  Petersen,  1974). 
However,  based  on  our  earlier  work,  we  can  currently  compare  the  three 


methods  for  a  vertical  point-force  source.  This  paper  makes  this 
comparison.  For  the  case  where  the  source  medium  is  Yucca  Flat  tuff  and 
the  propagation  medium  is  CIT109  and  including  a  vertical  component  LP 
LRSM  instrument,  we  find  that  the  CLEF  method  most  closely  matches  the 
maximum  peak-to-peak  amplitude  of  the  signal  generated  by  the  RT  method. 
However,  the  phase  distortion  inherent  in  the  CLEF  method  reduces  the 
period  associated  with  the  maximum  amplitude. 


REPRESENTATION  THEOREM  VALIDATIONS. 

The  theoretical  developement  of  the  Representation  Theorem  method 

for  calculating  fundamental  mode  Rayleigh  waves  from  complex  sources  has 

been  developed  in  a  series  of  papers  and  presentations  (e.g.  Harkrider 

et  al.,  1979;  Apsel  et  al.,  1980;  Harkrider,  1981;  and  Glover  and 

Harkrider,  1981a,  b).  Basically  the  method  consists  of  surrounding  the 

source  by  a  closed  surface  V  over  which  the  stresses  and 

SES 

displacements  generated  by  the  interaction  of  the  source  with  the  source 
region  medium  can  be  monitored.  These  quantities,  which  we  refer  to  as 
forcing  functions,  are  then  convolved  with  the  appropriate  Green's 
functions  for  the  exterior,  or  propagation  medium.  Details  of  this 
proceedure  and  corresponding  equations  can  be  found  in  the  references. 
As  we  have  previously  pointed  out,  the  RT  method  permits  the  convenient 
extension  of  detailed  calculations  of  complex  source/structure 
interactions  to  distances  at  which  meaningful  estimates  of  Mg  can  be 
made. 

Harkrider  (1981)  presented  some  preliminary  results  in  which  we 
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compared  the  surface  wave  generated  by  a  buried  vertical  point-force  in 
a  homogenous  half space  evaluated  using  the  RT  method  with  that 
calculated  from  the  analytic  formulation  of  Karkrider  (1964,  1970).  The 
results  in  both  cases  were  somewhat  "saw-toothed"  due  to  problems 
associated  with  picking  the  origin  time  of  the  synthetic  time  series. 
We  have  since  improved  our  numerical  techniques  for  evaluating  the 
Green's  functions  and  eliminated  this  problem.  Before  we  evaluate  the 
mixed  path  techniques,  we  will  show  these  new  results  in  order  to 
establish  the  accuracy  of  the  RT  method. 

Figure  1  shows  the  geometry  for  the  homogenous  halfspace  run.  £ 

SES 

is  a  cylinder  of  revolution  with  height  and  radius  2.1km.  There  are  11 

nodes  per  side  with  a  spacing  of  0.2km  (the  corner  node's  contribution 

is  evaluated  twice).  The  source  is  placed  at  0.4km  depth.  The  forcing 

functions  are  then  evaluated  at  each  node  using  the  finite  element  code 

SWIS  (Frazier  and  Petersen,  1974).  Figure  2  shows  the  contribution  to 

the  final  solution  for  the  sum  of  the  nodes  on  the  bottom  and  the  side 

of  T  ,  identified  by  forcing  function  component,  whereas  Figure  3 
SES 

shows  the  corresponding  sums  for  all  the  nodes.  The  total  RT  result  is 
also  shown  in  Figure  3  together  with  the  analytic  (Direct)  result.  All 
time  series  are  sampled  at  5  samples/sec  and  a  cosine  taper  has  been 
applied  over  the  entire  spectrum.  The  amplitude  of  the  RT  result  is  91Z 
of  the  Direct  result,  and  the  contributions  from  the  stresses  and 
displacements  are  approximately  equal.  However,  when  the  results  are 
convolved  with  a  vertical  component  LP  LRSM  seismometer  response 
(Figures  4  and  S)  the  final  RT  solution  is  dominated  by  the  contribution 
from  the  stresses,  specifically  the  normal  stress  on  the  bottom  and  the 
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tangential  stress  on  the  side  in  the  approximate  ratio  2:1. 

The  material  properties  of  the  homogenous  half space  run  (Figure  1) 
are  almost  Identical  to  those  of  the  upper  14km  of  the  layered  model 
CIT109  (Table  1).  Therefore  we  can  use  the  forcing  functions  from  the 
former  with  the  Green's  functions  from  the  latter  without  introducing 
any  appreciable  error.  Figures  6  and  7  show  the  corresponding  results, 
again  for  a  LP  LRSM  instrument,  sampled  at  5  samples /sec,  but  at  1200  km 
epicentral  distance  and  without  any  tapering.  In  this  case  the 
peak-to-peak  amplitude  of  the  RT  result  is  972  of  the  amplitude  of  the 
direct  result  and  the  associated  period  of  13.6  sec  is  identical  in  both 
cases.  For  the  dispersed  signal,  the  final  result  is  clearly  dominated 
by  the  contribution  from  the  vertical  component  of  stress  acting  on  the 
bottom  and  side  of  the  cannlster. 

ANALYTIC  APPROXIMATIONS. 

The  first  of  the  analytic  approximations  that  we  compare  using  the 
RT  method  is  the  CLEF  method  applied  by  Bache  et  al.  (1978)  in  an 
analysis  of  NTS  events  recorded  at  Tuscon,  Arizona  and  Albuquerque,  New 
Mexico.  In  Harkrider's  notation,  the  frequency  domain  vertical 
component  response  to  a  point  force  in  the  mixed  media  may  be  written  as 


The  second  method  is  that  used  by  Alexander  et  al.  (1976)  in  an 


analysis  of  a  sequence  of  earthquakes  at  Oroville,  California  recorded 
at  the  SRO  station  at  Albuquerque.  The  vertical  component  for  this 
approximation  is  given  by 


where  the  subscript  s  and  r  denote  properties  of  the  source  and 
propagation  median  respectively.  These  expressions  can  be  evaluated 
using  slight  modifications  of  computer  codes  used  to  calculate  the 
Direct  results  shown  in  Figures  5  and  7  with  little  appreciable  increase 
in  cost. 

MIXED  PATH  CALCULATIONS. 

With  a  minor  modification  of  the  finite  element  mesh,  we  can 
examine  the  Rayleigh  waves  from  a  point-force  in  a  medium  that  differs 
from  that  of  the  receiver.  Figure  8  shows  the  geometry  for  the  runs 
that  we  made  where  the  receiver  was  in  the  same  homogenous  propagation 
medium  as  in  Figure  1.  The  source  region  medium  was  a  cylinder  with  a 
height  and  radius  of  1.8km  with  material  properties  of  either  Yucca  Flat 
tuff  or  Climax  Stock  granite  (Table  1).  By  keeping  the  position  of  Y 

SES 

fixed,  we  were  able  to  monitor  the  first  6  sec  of  energy  transmitted 
across  the  Impedance  boundary. 

Figure  9  shows  the  results  for  the  two  source  media  typical  of  NTS. 
For  the  Climax  Stock  case  the  rigidity  ratio  is  approximately  0.6  and 


the  waveforms  given  by  the  three  methods  are  similar.  However  the 
difference  in  amplitudes  is  significant.  For  the  Yucca  Flat  case,  where 
the  rigidity  ratio  is  approximately  0.1,  the  waveforms  from  the  three 
methods  are  quite  dissimilar  and  the  amplitude  variation  is  large.  Both 
analytic  results  broaden  the  signal  width  compared  to  the  RT  result, 
indicating  some  phase  distortion.  When  the  traces  are  convolved  with 
the  response  of  a  LP  LRSM  instrument  (Figure  10),  the  differences 
between  results  of  the  three  computational  methods  are  noticeably 
reduced.  For  Climax  Stock,  the  amplitudes  given  by  the  CLEF  and  UTC 
method  are  identical,  though  still  somewhat  larger  than  the  RT  result. 
For  Yucca  Flat  ,  the  amplitude  disparities  between  the  three  methods  are 
considerably  less.  The  CLEF  result  shows  some  differences  with  respect 
to  the  RT  result  at  the  higher  frequencies,  but  not  so  great  as  the  UTC 
result  which  is  dominated  by  high  frequency  energy. 

We  also  convolved  the  Climax  Stock  forcing  functions  with  the 
CIT109  Green's  functions  for  a  vertical  component  LP  LRSM  intsrument  at 
1200km.  Figure  11  shows  the  RT  method  result  together  with  the  analytic 
results.  The  peak-to-peak  amplitudes  of  the  RT  and  CLEF  results  closely 
agree.  The  associated  periods  are  15  sec  and  12  sec  respectively.  The 
UTC  method  gives  amplitudes  a  factor  of  3  too  large  at  periods  around  6 
sec  and  these  distort  the  signal.  The  amplitude  of  the  13  sec  period  is 
approximately  equal  to  that  given  by  the  other  methods. 

It  is  interesting  to  note  that  the  peak-to-peak  amplitude  of  the 
mixed  path  result  computed  using  the  RT  method  is  4.0x10“^*  cm  compared 
to  3.6x10”^  cm  for  the  pure  CIT109  result.  The  source  in  both 
Instances  had  a  strength  of  1  dyne. 
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DISCUSSION. 

There  are  two  main  conclusions  that  can  be  drawn  from  these 
results.  These  are: 

i)  In  order  to  compute  the  RT  result  it  is  only  necessary  to 
perform  the  convolutions  using  the  vertical  component  of  stress  on  £ 

SES 

when  the  result  is  to  be  put  through  a  LP  LRSM  instrument  in  order  to 
compute  an  Mg.  This  should  effect  a  considerable  saving  in 
computational  effort  when  evaluating  both  forcing  and  Green's  functions. 

li)  that  for  the  type  of  mixed  paths  considered  here,  and  when  seen 
through  a  LP  LRSM  instrument,  the  CLEF  method  gives  results  comparable 
to  the  RT  method. 

In  light  of  these  results,  we  are  undertaking  a  similar  analysis 
using  a  finite  element  simulation  of  an  explosion  source  in  a  half  space 
to  generate  forcing  functions. 
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Figure  Captions 

Run  geometry  for  a  vertical  point  force(s)  in  a 

homogeneous  half-space.  X's  denote  nodes  on  £  at 

SES 

which  forcing  functions  are  calculated. 

Contributions  from  sum  of  modes  on  bottom  and  side  of 

Y  ,  Identified  by  forcing  function  component,  for 
SES 

problem  shown  in  Figure  1.  Amplitudes  are  ground 

displacement  in  cm  x  10“^  due  to  force  of  1  dyne. 

Total  contributions,  R.  T.  and  Direct  solutions. 

Units  as  in  Figure  2. 

Results  as  in  Figure  2,  but  convolved  with  LPZ  LRSM 
response,  plotted  on  common  scale. 

Results  of  Figure  3  after  convolution  with  LPZ  LRSM 
instrument  response. 

Contributions  from  bottom  and  side  of  7  for  vertical 

SES 

point  force  of  1  dyne  in  model  CIT109,  convolved  with 
LPZ  LRSM  instrument  of  1200  km  distance.  Units  are  cm 
x  10"15. 

Total  contribution,  RT  and  Direct  solutions  for  model 
CIT109. 

Geometry  for  mixed-path  runs  shown  in  Figures  9  and  10. 


For  Figure  11,  receiver  R  was  placed  at  A  -  1200  km 


Figure  9.  Comparison  of  the  mixed-path  results  for  2  NTS  source 
media,  scale  factors  shown  separately. 

Figure  10.  Results  of  Figure  9  convolved  with  LPZ  LRSM  response. 
Figure  11.  Mixed-path  results  for  model  C1T109,  LPZ  LRSM 
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ABSTRACT 

Vaveform  modeling  of  teleseismic  long-period  body  wave  phases  of 
shallow  earthquakes  has  proven  quite  effective  in  determining  source 
parameters  for  events  larger  than  magnitude  six.  Unfortunately,  these 
teleseismic  phases  become  too  weak  for  the  smaller  events  and  regional 
data  mudt  be  used  which  are  generally  much  more  complicated.  This  is 
because  the  crust-mantle  system  is  acting  like  a  leaky  waveguide  and  a 
large  number  of  rays  are  required  to  model  the  observations.  In  this 
report  we  review  a  procedure  for  the  systematic  determination  of  source 
parameters  from  the  regional  body  waves.  A  least-squares  inversion 
technique  which  is  based  on  a  cross-correlation  of  the  data  and  a 
synthetic  is  used.  The  synthetics  are  a  linear  combination  of  the  three 
fundamental  faults.  A  set  of  profiles  of  the  synthetics  is  presented. 
Although  the  synthetics  are  generated  for  a  model  developed  for  the 
western  U.S.,  this  model  seems  to  be  adequate  for  most  continental 
regions.  The  inversion  is  parameterized  in  terms  of  fault  strike,  dip 
and  rake;  these  parameters  are  relatively  Insensitive  to  small  changes 
in  crustal  structure.  Several  examples  are  presented. 
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INTRODUCTION 

A  considerable  amount  of  effort  has  been  expended  to  determine  the 
source  parameters  of  moderate  size  earthquakes,  although  it  can  be 
besieged  with  difficulties.  Ideally,  a  large  amount  of  information  can 
be  derived  from  modeling  the  long-period  body  waves  (see  Helmberger, 
1974;  Langston  and  Helmberger,  1975).  Unfortunately,  if  the  earthquake 
is  too  small  to  be  well  recorded  teleseismically,  which  is  the  case  for 
many  events  in  the  magnitude  range  between  5  and  6,  the  fault  plane 
orientation  must  be  constrained  by  local  short-period  data  and  the 
seismic  moment  usually  can  not  be  determined  unambiguously.  The 
World-Wide  Standard  Seismograph  Network  (WWSSN)  supplemented  by  other 
long-period  stations  and  arrays  provides  sufficiently  dense  coverage  in 
that  most  moderate  size  earthquakes  occurring  in  continental  regions 
will  produce  some  on-scale  records  of  the  long-period  body  waves  at 
regional  distances.  In  the  regional  distance  range  (1  -12°)  the  wave 
guide  properties  of  the  crust  produce  complicated  body  wave  signals; 
however,  in  most  cases  the  long-period  waveform  is  quite  distinctive  and 
sufficiently  Insensitive  to  crustal  structure  details  to  allow  the 
separation  of  the  source  and  structural  information. 

In  this  paper  we  review  a  procedure  for  extracting  the  source 
parameters  of  moderate  size  earthquakes  from  long-period  regional 
phases.  The  technique  Involves  an  Iterative  inversion  process  which 
minimizes  the  difference  between  a  synthetic  seismogram  and  the 
observation.  The  synthetics  are  constructed  using  Green's  functions 
computed  for  a  single,  very  simple  structure.  These  Green's  functions 
appear  to  be  an  adequate  model  for  most  continental  regions,  thus 
allowing  a  quick  and  approximate  determination  of  the  fault  parameters. 
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The  inversion  is  parameterized  in  terms  of  strike  dip  and  rake.  The 
number  of  inversion  parameters  has  been  minimized  so  that  inadequacies 
in  the  Green's  functions  are  not  over  emphasised.  Obviously,  the 
structural  model  is  more  approplate  for  certain  region  than  others,  so 
the  inversion  parameters  chosen  are  those  which  are  most  robust.  The 
main  advantage  of  this  technique  is  that  it  only  requires  a  small  data 
set.  The  general  usefulness  of  this  technique  is  illustrated  by 
inverting  regional  data  from  earthquakes  occurring  in  the  United  States, 
northern  Canada  and  southern  Europe. 

THE  GREEN'S  FUNCTIONS 

The  techniques  for  constructing  the  Green's  functions  are  discussed 
in  detail  elsewhere  (Helmberger  and  Engen,  1980;  Wallace  et  al.,  1981). 
What  is  discussed  here  is  some  simplifying  approximations  and  their 
applicability.  The  waveform  of  interest  is  that  part  of  a  regional 
seismogram  which  arrives  before  the  S-wave.  This  waveform  is  referred 
to  as  Vyii*  It  Is  simplest  to  discuss  P^  in  terms  of  rays  which  are  in 
a  waveguide.  The  first  part  of  P^  is  dominated  by  P-waves  (Pn)  and, 
moving  back  into  the  record,  the  waveform  contains  progessively  more  SV 
(PL)  contributions.  The  SV  energy  corresponds  to  rays  which  are 
reflected  within  in  the  crust  and  undergo  subsequent  mode  changes  at  the 
free  surface  and  the  Moho.  The  interference  of  all  the  rays  gives  rise 
to  the  waveform;  parameters  such  as  crustal  thickness  or  velocity 
contrast  between  the  crust  and  Moho  control  the  waveform  dispersion. 
Figure  1  shows  some  typical  waveforms. 

The  Green's  functions  are  constructed  by  summing  generalized  rays 
for  a  point  shear  dislocation.  As  an  example,  consider  the  following 
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equation  for  the  vertical  displacement  in  cylindrical  coordinates: 

/ .  3  v 

w(r,z,6,t)  -  ^D(t)  *  £  W1(t)  A1J  (1) 

where  D(t)  is  the  far  field  time  history,  the  source  region  density 
and  M0  is  the  seismic  moment.  The  summation  process  adds  the 
contributions  of  the  three  fundamental  faults;  the  are  the  Green's 
functions  for  vertical  strike-slip,  vertical  dip-slip  and  45°  dip-slip 
step  dislocations  respectively.  The  A^  are  coefficients  determined  by 
source  orientation  and  are  given  by: 

(0,X,6)  -  sin  2G  cos  X  sin  6  +  1/2  cos  20  sin  X  sin£6 
A2  (0,X,6)  *  cos  0  cos  X  cos  6  -  sin  0  sin  X  cos  26  ^2) 

A3  (0, X, 6)  -  (1/2)  sin  X  sin  26 

where  8  is  the  receiver  azimuth  from  the  end  of  the  fault  plane,  X  is 
the  rake  angle  and  6  is  the  dip  angle.  The  total  displacement  is  then 
the  sum  of  the  displacements  from  each  ray.  The  number  of  rays  which 
are  used  is  determined  by  the  structure  and  the  pass  band  of  the 
observation.  In  our  previous  work  we  have  shown  that  a  single  layer, 
corresponding  to  the  crust,  over  a  halfspace  mantle  is  a  sufficient 
structural  model  for  regional  long-period  records  to  allow  the 
extraction  of  source  parameters  of  moderate  size  earthquakes. 

As  an  example,  Figure  2  shows  a  comparison  of  the  synthetics  and 
the  records  of  the  1966  Truckee,  California  earthquake  (which  will  be 
discussed  in  more  detail).  The  synthetics  were  constructed  from  Green's 
functions  computed  for  the  crustal  model  in  Table  1  and  the  fault 
orientation  was  determined  by  the  Inversion  of  the  regional  data.  In 
this  case  the  source  time  function  and  moment  were  determined  by  other 
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methods  so  that  the  synthetic  amplitudes  can  be  viewed  as  predictions. 
The  numbers  on  the  traces  are  maximum  peak  to  peak  amplitudes.  The  only 
noticeable  difference  between  the  data  and  the  synthetics  is  the  high 
frequency  content  which  could  be  caused  by  several  things;  (1)  the 
effects  of  attenuation  within  the  crust  have  not  been  added  to  the 
synthetics,  and  (2)  the  very  sharp  boundaries  in  our  model  are  efficient 
in  trapping  short-period  energy.  Overall,  the  fit  of  the  synthetics  to 
the  data  justifies  the  use  of  the  simple  model,  and  the  high  frequency 
content  of  the  synthetics  does  not  effect  our  ability  to  determine  the 
source  parameters. 

The  only  expensive  or  complicated  part  of  this  modeling  process  is 
the  actual  generation  of  the  Green's  functions.  For  this  reason. 
Figures  3  and  4,  which  give  profiles  of  the  vertical  and  radial 
responses  for  the  three  fundamental  faults  are  presented.  Most 
earthquakes  which  produce  an  on-scale  at  long-period  WWSSN  stations 
have  similar  time  functions  (which  is  simply  a  reflection  of  the  event 
size).  In  Figures  3  and  4  a  trapezoid  with  a  1  second  rise,  1  second 
top  and  1  second  fall  is  used  for  the  far-fleld  time  history.  The 
displacements  have  also  been  convolved  with  a  15-100  Instrument. 
Because  of  the  differences  in  high  frequency  content  between  the  data 
and  synthetics  the  displacement  responses  have  been  filtered.  The 
filter  has  an  Impulse  response  of  a  triangle  which  has  a  2  second  rise 
and  fall.  When  comparing  these  displacements  with  data,  the 
observations  should  be  similarly  filtered.  Once  the  response  of  the 
three  fundamental  faults  is  known  any  seismogram  can  be  constructed  by  a 
linear  combination  of  them. 

The  displacements  in  figure  3  and  4  were  computed  for  a  source 


depth  of  8  km.  Varying  the  source  depth  between  5  and  15  km  has  only  a 
small  effect  on  the  waveform.  This  is  easily  understood  by  considering 
that  to  first  order,  a  change  in  source  depth  only  affects  the  travel 
time  of  the  first  segment  of  any  ray.  Figure  5  shows  a  comparison  of 
the  synthetics  at  1000km  for  three  different  source  depths.  After 
doubling  the  source  depth  (from  8  to  16  km)  the  essential  character  of 
the  waveform  is  still  preserved  and  the  source  information  is 
retrlvable.  In  contrast,  a  similar  change  in  crustal  thickness  would 
affect  the  travel  time  of  each  leg  of  a  given  ray,  hence  significantly 
changing  the  waveform  dispersion  (Wallace  and  Helmberger,  1980).  The 
insensitivity  of  the  displacements  to  source  depth  allows  the  responses 
in  figures  3  and  4  to  be  used,  at  least  in  a  qualitative  fashion,  to 
determine  the  source  parameters  of  most  crustal  earthquakes. 

The  only  other  major  question  of  applicability  of  the  presented 
displacements  is  the  structure  used  in  their  calculation.  Experience 
indicates  that  the  simplistic  model  is  justified.  The  model  in  Table  1 
is  an  average  developed  for  the  western  U.S.  although  it  appears  to  be 
sufficient  for  most  continental  regions  in  the  world.  The  waveform 
dispersion  is  dependent  on  crustal  thickness  and  the  contrast  between 
the  upper  mantle  P-velocity  and  the  mean  crustal  P-velocity,  so 
obviously  in  regions  with  anomalous  crustal  structure  such  as  the 
Tibetan  Plateau  the  responses  in  figures  3  and  4  would  be  inadequate. 
Also,  the  use  of  the  half space  to  approximate  the  upper  mantle  must 
break  down  at  some  point;  at  some  distance  a  significant  amount  of 
energy  will  be  present  in  the  form  of  diving  rays  which  have  turned  in 
the  mantle.  For  most  continental  regions  this  distance  appears  to  be 
about  12°.  If  diving  rays  are  present  the  ratio  of  the  Pn  to  PL 


amplitude  should  differ  from  the  synthetics  (Wallace  and  Helmberger, 
1981).  Also,  the  predicted  amplitude  should  begin  to  diverge  from  the 
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observed  due  to  attenuation  in  the  upper  mantle.  Figure  6  shows  a 
profile  of  the  Truckee  records;  the  amplitudes  have  been  corrected  for 
the  azimuthal  radiation  pattern  and  the  polarln'-'  have  been  adjusted  to 
show  a  smooth  dispersion  pattern.  To  the  right  is  a  profile  of 
displacements  for  a  strike-slip  earthquake  (such  a  profile  can  be 
constructed  from  figure  3).  Note  that  there  does  not  appear  to  a 
systematic  break  down  in  waveform  shape  or  amplitude  over  the  distance 
range  to  4°  to  12°.  Also  note  that  stations  BOZ  and  TUC  are  nearly 
nodal  and  thus  their  amplitudes  are  not  particularly  reliable. 

INVERSION  TECHNIQUE 

The  ability  to  determine  the  source  parameters  of  an  earthquake  by 
comparing  an  observed  with  a  predicted  waveform  depends  on  the 
assessment  of  the  quality  of  fit.  In  a  previous  paper  (Wallace  et  al., 
1981)  we  presented  a  least-squares  waveform  inversion  technique  which 
makes  use  of  an  error  function  determined  by  the  cross-correlation  of  a 
long-period  seismogram  and  a  synthetic; 


where  f  is  the  observed  record,  g  is  the  synthetic  and  the  Integral  is  a 
zero  lag  cross-correlation.  The  limits  of  integration  are  the  time 
length  of  the  window  in  which  the  waveforms  are  correlated.  The 
denominator  serves  to  normalize  both  the  data  and  the  synthetics.  This 
normalization  process  makes  the  error  function  insensitive  to  the 
absolute  amplitudes.  To  minimize  the  error,  which  corresponds  to 
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maximizing  the  correlation,  we  allow  f  and  g  to  optimally  align 
themselves  with  regard  to  waveform.  f  and  g  are  aligned  a  priori  in 
time  by  matching  first  breaks  and  ignoring  absolute  travel  time.  The 
error  function  can  be  rewritten  by  considering  that  the  synthetic 


seismogram 

can  be 

constructed  with  the 

three 

fundamental 

faults. 

In 

this  case 

there 

is 

a 

summation  of 

cross-correlations 

between 

the 

observed  and  each 

of 

the 

fundamental  faults. 

For  a  given 

range , 

the 

croDS-correlations  are  constant  and  errors  can  be  minimized  by  varying 
the  A^'s  in  equation  2.  In  other  words,  once  the  cross-correlations  are 
computed  the  source  orientation  is  determined  iteratively  and  only  the 
constants  have  to  be  recalculated.  For  a  detailed  discussion  of  the 
inversion  technique  see  Vallace  et  al.  (1981). 

Once  the  source  orientation  is  fixed,  the  moment  of  an  earthquake 
can  be  determined  by  comparing  the  amplitude  of  the  synthetics  and 
observations.  Adopting  the  units  of  Helmberger  and  Malone  (1975),  and 
expressing  the  range  in  km,  time  in  seconds,  density  in  gm/cm^,  velocity 
in  km/sec,  the  moment  in  dyne-cm  and  displacement  in  cm  yields: 

M  -  4xo  x  1020  (  data  amplitude  .  (4) 
o  '  synthetic  amplitude 

A  moment  can  be  determined  by  comparing  the  maximum  peak  to  peak 
amplitude  for  any  time  window  used  for  the  correlation.  It  has  been 
found  that  a  moment  should  be  determined  for  a  few  peaks  at  a  given 
station.  The  ratio  of  the  moment  at  each  station  to  the  mean  is  a 
measure  of  the  amplitude  stability.  In  general,  the  moments  determined 
from  are  in  very  good  agreement  with  those  determined 
teleselsmically,  using  an  assumption  of  t*  ■  1  sec. 
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EXAMPLES 

We  have  inverted  the  waveform  from  five  earthquakes  to 

demonstrate  the  utility  of  the  technique.  Three  of  the  earthquakes  are 
in  the  western  U.S.  Two  other  earthquakes,  one  in  Baffin  Bay  in  the 
Artie  and  the  other  in  Turkey,  have  been  included  to  demonstrate  that 
the  Green's  functions  are  not  unique  to  the  U.S.  Both  dip-slip  and 
strike-slip  mechanisms  are  represented  in  the  suite  of  examples. 

Truckee,  California  (9/12/67) 

The  Truckee  earthquake  was  a  strike-slip  event  at  10  km  depth  which 
produced  excellent  regional  records  but  very  few  teleseismlc  body  wave 
records  as  typical  of  moderate  size  strike-slip  events.  The  Truckee 
earthquake  (mb«5.7)  has  been  studied  by  numerous  authors  (Ryall  et  al., 
1968;  Tsai  and  Aki,  1970;  Burdick,  1977)  making  it  a  good  test  case. 
Tsai  and  Aki  (1970),  from  first  motion  studies  and  modeling  of  the 
surface  waves,  determined  this  event  to  be  pure  strike-slip  on  a  fault 
plane  striking  N44°E  and  dipping  80°SE.  The  surface  wave  moment  was 
determined  to  be  0.83  x  102^  dyne-cm.  Figure  2  shows  the  location  of 
the  epicenter,  recording  stations  and  filtered  data  for  Truckee.  Also 
shown  are  the  synthetics  determined  from  the  inversion  results.  Note 
that  BOZ  and  TUC  are  very  nearly  nodal.  The  inversion  yields  a 
mechanism  which  is  very  similar  to  Tsai  and  Aki's  (1970);  a  strike  of 

N43°E,  a  dip  of  76°SE  and  a  rake  of  -11°.  The  only  significant 

difference  is  the  slight  dip-slip  component  in  our  solution,  which  is 
also  acceptable  on  the  basis  of  the  first  notion  data.  The  moment 

determined  from  the  waveforms  is  0.87  x  1025  dyne-cm,  which  is  in 


excellent  agreement  with  Tsai  and  Aki  (1970). 


El  Golfo,  Mexico  (8/7/66) 

The  El  Golfo  earthquake  (mb“6.3,  Ms*6.3)  is  a  strike-slip  event 
which  occurred  near  the  mouth  of  the  Colorado  River  at  the  northern  end 
of  the  Gulf  of  California.  Ebel  et  al.  (1978)  determined  the  fault 
plane  to  be  striking  E140°S,  dipping  85°  to  the  southwest  and  with  a 
rake  of  183°,  and  determined  the  depth  of  the  event  to  be  10  km.  Using 
teleselsmic  long-period  P-waves  they  determined  a  moment  of  5.0  x  102-* 
dyne-cm. 

El  Golfo  is  about  the  maximum  size  event  which  can  be  used  in  the 
Inversion  technique.  The  P ^  records  are  barely  on  scale  at  the 
stations  used.  Figure  7  shows  the  location  of  the  epicenter,  recording 
stations  and  the  waveforms.  In  this  case  the  time  function  is  a 
triangle  with  a  two  second  rise  and  fall;  the  long-period  impluse  is  a 
reflection  of  the  event  size.  The  long-period  time  function  allows  us 
to  dispense  with  using  a  filter.  In  figure  7,  shown  below  the  observed 
waveforms  are  the  synthetics  for  the  inversion  solution;  a  strike  of 
E137°S,  a  dip  of  81°  and  a  rake  of  175°.  The  inversion  solution  is  in 
good  agreement  with  Ebel  et  al.  (1978).  The  moment  determined  from  the 
P  is  4.6  x  1025  dyne-cm. 

Orovllle,  California  (8/7/75) 

The  Orovllle  earthquake  (M«5.6)  was  a  normal  faulting  event  and  is 
interesting  because  the  surface  wave  (Hart  et  al. ,  1977)  and  body  wave 
(Langston  and  Butler,  1976)  analysis  yield  substantially  different 
moments.  Langston  and  Butler  (1976)  determined  the  strike  to  be  180°, 
with  a  dip  of  65°  and  a  rake  of  -70°.  Their  moment  determination  is  5.7 
x  102*  dyne-cm.  Hart  et  al.  (1977)  suggests  that  the  surface  waves  are 
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consistent  with  the  body  wave  mechanism  but  the  moment  is  a  factor  of  3 
larger  (1.9  x  10  dyne-cm).  Figure  8  shows  the  location  of  the  event, 
the  stations  used  in  the  inversion  analysis  and  the  filtered  data  and 
synthetics.  The  inversion  solution  has  shifted  the  mechanism  to  a 
strike  of  204°,  dipping  66°  and  a  rake  of  -85°.  This  new  solution  only 
violates  a  few  first  motions,  but  the  aftershock  trend  tends  to  support 
the  180°  strike.  The  moment  determined  from  is  6.9  x  10^  dyne-cm. 

Baffin  Bay,  Canada  (9/4/63) 

The  Baffin  Bay  earthquake  (M«5.9)  is  a  normal  event  associated  with 
a  continental  margin  area.  The  travel  path  to  each  of  the  stations  used 
in  the  Inversion  includes  portions  of  crustal  and  oceanic  regions  which 
makes  an  ideal  event  to  test  the  applicability  of  the  Green's  functions. 
Liu  and  Kanamori  (1980)  modeled  the  body  waves  and  determined  a  fault 
plane  solution  with  a  strike  of  98°,  dip  of  66°N  and  a  rake  of  -103°. 
The  location  of  the  event  and  the  filtered  data  and  synthetics  are  shown 
in  Figure  9.  The  Inversion  solution  has  a  mechanism  striking  74°, 
dipping  66°  and  has  a  rake  of  -100°.  Again  the  only  parameter  which  is 
appreciably  different  from  the  teleseismic  analysis  is  the  strike.  In 
any  case  the  Inversion  solution  is  acceptable,  and  the  difference  of  20° 
in  strike  can  be  considered  the  resolution  for  dip-slip  events. 

Turkey  (6/13/65) 

The  Turkish  event  (M*5.1)  is  a  shallow  normal  event  which  occurred 
In  southwest  Turkey  in  a  region  of  north-south  extension.  McKenzie 
(1972)  used  first  motion  data  to  determine  a  pure  normal  mechanism  with 
a  strike  of  101°,  dipping  70°  to  the  south,  although  it  is  not  well 
constrained.  There  were  three  WWSSN  stations  at  regional  distances 
which  could  be  used  in  the  Inversion  process.  Figure  10  shows  the 
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location  of  the  event  and  the  recording  stations.  The  filtered  data  and 
the  fit  of  the  synthetics  are  also  shown.  The  inversion  solution 
(strike  ■  131°,  dip  *  68°,  rake  ■  -88°)  is  conrlstent  with  the  first 
motion  data,  although  it  differs  in  strike  from  McKenzie's  solution. 
Again  the  three  station  solution  is  quite  acceptable,  considering  the 
quality  of  the  first  motion  data. 


DISCUSSION 

Determining  the  fault-plane  orientation  of  moderate  size 
earthquakes  is  often  a  frustrating  experience  due  to  the  paucity  of  high 
quality  data.  Earthquakes  in  the  magnitude  range  of  5  to  6  are  quite 
Important  and  often  are  the  only  "measurable"  expression  of  the  present 
tectonic  environment.  All  the  available  data  must  be  used  to  extract 
the  source  parameters  of  these  moderate  size  events,  and  the  modeling  of 
waveforms  can  provide  a  valuable  constraint  in  this  process.  Every 
situation  will  probably  be  unique  and  it  is  difficult  to  predict  which 
data  set  will  be  the  most  definitive.  Nevertheless,  it  appears 
worthwhile  to  consider  the  inversion  seperately  and  its 

resolvability  dependence  on  source  orieutation. 

The  lnverison  process  we  discussed  only  requires  a  small  data  set. 
With  ideal  azimuthal  separation,  a  data  set  comprised  of  Just  3  stations 
(vertical  and  radial  component)  can  yield  good  solutions  (it  is  possible 
to  construct  a  case  where  the  inversion  is  unstable,  but  in  practice 
this  has  never  happened).  In  almost  all  cases  4  recording  stations  are 
sufficient.  ^  yJL  should  be  polarized  in  the  vertical  and  radial  planes. 
Rotation  of  the  horizontal  components  for  a  number  of  events  Indicates 
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that  there  is  very  little  energy  on  the  tangential  component  implying 
little  contamination  from  such  effects  as  multipathing.  Therefore,  it 
is  usually  sufficient  to  take  the  largest  of  the  horizontal  components 
to  be  the  radial  waveform  in  the  inversion.  The  resolving  ability  of 
the  Inversion  (or  conversely  the  error)  depends  on  the  type  of 
earthquake.  The  experience  gained  by  considering  the  examples  presented 
in  the  last  section  indicates  that  the  mechanisms  of  strike-slip 
earthquakes  can  be  determined  quite  well  with  relatively  few  stations. 
The  strike  is  usually  determined  to  within  5°  of  that  determined  by 
other  methods.  The  rake  is  the  least  resolvable  parameter  for 
strike-slip  events  and  can  vary  up  to  15°  from  that  determined  by  first 
motion  studies.  The  mechanisms  of  dip-slip  earthquakes  are  more 
difficult  to  determine.  Although  the  dip  and  rake  are  usually 
determined  in  good  agreement  with  other  studies,  the  strike  may  vary  up 
to  20°.  This  feature  is  illustrated  by  considering  a  45°  dipping  normal 
fault.  In  this  case,  most  regional  stations  lie  within  the 
compressional  region  of  the  focal  sphere  and  any  given  azimuth  will 
produce  remarkably  similar  waveforms.  Fortunately,  dip-slip  events  have 
rather  strong  teleselsmic  P-waves  but  again  the  waveforms  are  all  the 
same  with  little  dependence  on  azimuth.  In  this  case  the  stations  lie 
in  the  center  of  the  focal  sphere  and  are  all  dllatational.  As  an 
example  consider  Figure  11,  which  displays  the  focal  sphere  and  the 
teleselsmic  waveforms  for  the  Oroville  earthquake.  Note  that  the 
P-waveforms  are  all  similar  and  the  synthetics  (Langston  and  Butler, 
1976)  do  not  add  much  insight  into  determining  the  strike  of  the  fault. 
Figure  12  shows  the  filtered  regional  data  and  synthetics  computed  from 
the  teleselsmic  fault  parameter  determinations.  A  comparison  of  figure 
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12  and  8  will  show  that  the  regional  data  inversion  solution  improves 
the  fit  of  the  synthetics,  in  particular  PAS.  This  suggests  that  a 
logical  approach  would  be  to  invert  some  teleseismic  data  and  the 
regional  data  simultaneously.  Since  the  inversion  technique  relies  on 
the  cross-correlation  of  data  and  a  synthetic  the  joint  inversion  is 
quite  tractable. 

The  higher  resolution  of  the  strike-slip  events  is  actually 
fortuitous.  Moderate  size  strike-slip  earthquakes  rarely  produce  usable 
teleseismic  P-waves  due  to  their  inefficiency  in  radiating  energy 
straight  down.  On  the  other  hand,  strike-slip  events  produce  very  good 
regional  waveforms.  This  allows  the  inclusion  of  a  larger  data  set  in 
the  inversion  and  hence,  the  resolution  problem  is  at  least  partially 
resolved. 

It  is  reasonable  to  consider  what  effect  the  structual  model  has  on 
the  inversion  results.  As  a  test  of  the  insensitivity  of  the  fault 
orientation  to  small  changes  in  crustal  parameters  the  El  Golfo 
earthquake  was  reinverted  with  a  different  structure.  The  crustal 
thickness  was  reduced  to  24  km,  the  source  depth  was  moved  to  12  km  and 
the  Pn  velocity  reduced  to  7.8  km/sec.  Although  the  quality  of  the  fit 
decreases  significantly  the  mechanism  returned  by  the  inversion  is 
similar;  strike  is  138°,  dip  is  82°  and  a  rake  of  181°.  The  moment 
increases  to  6.9  x  10^  dyne-cm. 

CONCLUSIONS 

It  is  possible  to  extract  the  source  parameters  of  moderate  size 
earthquakes  from  the  long-period  regional  body  waves.  The  procedure 
requires  the  comparison  of  the  observed  waveform  with  a  synthetic;  the 
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synthetics  can  be  generated  by  a  linear  combination  of  the  waveforms  of 
the  three  fundamental  faults  shown  in  figures  3  and  4.  Although  these 
synthetics  are  for  a  single  model  the  inversion  parameters  (fault 
strike,  dip  and  rake)  are  fairly  insensitive  to  small  changes  in  crustal 
thickness,  Pn  velocity  and  mean  crustal  velocity.  This  allows  this 
single  set  of  Green's  functions  to  be  used  for  most  continental 
earthquakes.  The  Inversion  procedure  only  requires  a  small  data  set, 
and  is  particularly  ideal  for  strike-slip  earthquakes. 
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Figure  1:  The  horizontal  component  of  motion  for  four  moderate  size 
earthquakes.  The  magnitude  and  distance  to  the  recording  station  is 
given  for  each.  Clockwise  from  the  upper  left;  (1)  Sept.  12,  1966, 
Truckee,  California,  a  strike-slip  event,  (2)  June  13,  1965,  southwest 
Turkey,  a  normal  event,  (3)  Aug.  7,  1975,  Oroville,  California,  a 
normal  event,  and  (4)  Dec.  10,  1967,  off  the  coast  near  Cape  Mendocino, 
California,  a  strike-slip  event. 

Figure  2:  The  vertical  waveforms  of  the  1966  Truckee  earthquake. 
The  star  denotes  the  epicenter.  The  data  is  the  top  trace  at  each 
station  and  the  trace  below  is  the  synthetic  fit.  The  strike-slip 
mechanism  has  two  nodal  planes  which  project  through  TUC  and  BOZ.  To 
the  right  of  each  trace  is  the  observed  or  predicted  amplitude  (on  the 
basis  of  .8  x  102^  dyne-cm)  in  10”^  cm. 

Figure  3:  Theoretical  displacement  profiles  for  the  vertical  component. 
The  Green's  fuctions  were  computed  from  the  model  presented  in  Table  1 
and  have  been  convolved  with  a  source  time  function  represented  by  a 
trapezoid  (ti*»l,t2“l,t3"l),  a  triangular  filter  (2  second  rise  and 
fall),  and  a  WWSSN  long-period  instrument. 

Figure  4:  Displacement  profiles  for  the  radial  component.  The  Green's 
functions  are  computed  every  100  km.  They  have  been  convolved  with  the 
time  function,  Instrument  and  filter  discribed  in  figure  3. 

Figure  5:  V  waveforms  at  1000  km  for  three  different  source  depths. 

Figure  6:  The  Truckee  waveforms  have  been  corrected  for  horizontal 
radiation  pattern  and  plotted  as  a  function  of  distance.  The  maximum 
amplitude  is  shown  to  the  right  of  each  trace.  Mote  that  BOZ  and  TUC 
are  very  close  to  nodes. 

Figure  7:  Data  and  synthetics  from  the  El  Golfo  earthquake.  The  map 
gives  the  location  of  the  event  (6tar)  and  the  recording  stations. 
Along  each  trace  is  the  ratio  of  the  station  moment  and  the  average 
moment. 

Figure  8:  Filtered  data  and  synthetics  from  the  Oroville  earthquake. 
At  all  the  stations  except  G0L  both  the  vertical  (the  first  trace  pair) 
and  radial  components  are  shown. 

Figure  9:  Location  of  the  Baffin  Bay  earthquake  (star)  and  the 
recording  stations.  The  filtered  data  and  synthetics  from  both  the 
vertical  and  radial  components  are  shown. 

Figure  10:  Location  of  the  Turkey  event  (star)  and  the  recording 
stations.  The  filtered  data  and  synthetics  for  both  the  vertical  and 
radial  components  are  shown. 

Figure  11:  The  teleseismlc  waveforms  from  the  Oroville  earthquake. 
Note  the  similarity  of  the  vavefoms  at  all  azimuths.  Shown  on  each 
trace  is  the  moment  (x  1025  dyne-cm)  determined  at  that  station.  Note 


the  factor  of  3  scatter  (figure  from  Langston  and  Butler,  1976). 

Figure  12:  Filtered  P^from  the  Oroville  earthquake.  The  synthetics 
were  computed  with  the  teleselsmic  fault  plane  solution. 
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ABSTRACT 

A  numerical  method  for  evaluating  the  Klrchhoff-Helmholtz  Integral 
Is  described.  The  Klrchhoff  response  is  calculated  by  discretizing  the 
surface,  specifying  simple  point  sources  on  each  element  of  the  surface, 
and  summing  the  contribution  from  the  elements.  The  results  of  the 
method  are  compared  to  those  of  an  asymptotic,  first  motion 
approximation  of  the  analytical  solution  of  SH  waves  impinging  on  a 
rigid  sphere.  The  agreement  between  the  results  of  the  two  methods  is 
excellent  for  source  and  receiver  distances  which  are  large  compared  to 
the  radius  of  the  sphere.  The  method  is  applied  to  the  calculations  of 
reflections  from  mountain  topography  and  a  planar  surface  with  an 
aperture.  The  phase  shifts  of  pulses  are  consistent  with  optics;  the 
amplitudes  are  not.  The  method  does  predict  frequency  dependence  of  the 
scattered  amplitudes.  Calculations  are  presented  to  model  spall  which 
produce  travel-time  and  amplitude  anomalies  consistent  with  observations 


from  nuclear  blasts 


Introduction 


Many  wave  propagation  phenomena  cannot  be  adequately  modeled  by 
existing  solutions  to  plane-layered  media.  Tet  the  Increasing  use  of 
broad-band  seismic  data  to  determine  source  dislocations,  Q,  and 
velocity  structure  requires  a  knowledge  of  effects  of  material 
irregularities  in  the  medium  on  seismic  wave  propagation.  Certainly, 
documented  amplitude  and  dT/dA  anomalies  of  teleseismic  arrivals  at 
large  arrays  (Glover  and  Alexander,  1969;  Walck  and  Minster,  1980) 
which  vary  as  a  function  of  azimuth  suggest  the  existence  of  non-flat 
boundaries  at  depth. 

Numerical  schemes  which  handle  material  irregularities  are  in 
abundance.  Finite  difference  and  finite  element  codes  have  been  used 
successfully  (Boore,  1971;  Smith,  1975)  and  can  be  applied  to  a  variety 
of  materials;  however,  the  expense  of  calculating  the  response  at 
distances  which  are  large  compared  to  the  wavelength  of  interest  is 
prohibitive.  Rayleigh-FFT  techniques  have  been  exploited  for  these 
problems  (Aki  and  Lamer,  1970).  Similiarly,  Implementation  of  these 
methods  for  analysis  of  three  dimensional  scattering  is  also  costly. 
Geometric  ray  methods  are  useful  for  predicting  scattering  of  signals 
which  have  wavelengths  that  are  short  compared  to  the  size  of  the 
heterogeneity  (Hong  and  Helmberger,  1978).  But,  existing  ray  methods  do 
not  predict  frequency-dependent  amplitudes  of  scattered  pulses  and  do 
not  handle  diffracted  arrivals. 
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Discussion  of  techniques 

The  method  presented  In  this  paper  16  an  Integral  equation  approach 
and  is  based  on  the  evaluation  of  the  scalar  Integral  equation  variously 
called  the  Kirchhoff,  Helmholtz,  or  Huygen's  integral.  This  equation  is 
a  formulation  of  the  wave  equation  in  terms  of  a  linear  surface  Integral 
over  the  boundary  of  a  continuous  volume.  That  is,  the  scalar  wave 
equation  is 

2 

X  ~  <r,t)  -  V2u(r,t)  -  »(ro,t>  (1) 

at 


Here  u(r,t)  is  the  field  at  a  point  r  resulting  from  a  source  potential 
♦(rQ,£)  and  c  is  the  wave  speed.  Following  the  formalism  discussed  by 
Mow  and  Pao  (1971),  we  consider  the  motion  of  a  homogeneous  body  V  with 
a  smooth  boundary  3V  with  outward  pointing  normal  &.  then  if  £  e  V  and 
t  e  (0,-) 


dV 
o  o 


’  j(/0  '<£•  V  t-t.>»<£.-t>4t, 
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Here  G(£,£0,t-te)  is  the  fundamental  singular  solution  of  the  scalar 
wave  equation 


1  a  G  _2_. 
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Let  us  define  f(r^,t)  as  the  sum  of  the  first  two  Integrals  in  equation 
(2).  Then  f(£,t)  can  be  Interpreted  as  the  whole-space  solution  of  the 
problem  with  sources  and  initial  values  (rQ,0)  and  u(xo,0). 

Hence  u(x,t)  is  a  sum  of  the  direct  pulse  and  a  reflected  pulse  from  the 
surface  which  is  described  by  the  third  integral  in  equation  (2). 

If  jr  i  V,  then  the  left-hand  side  of  equation  (2)  is  aero.  If  X  e 
3V  then 

G(r,ro,t-t0)  J 

(4) 
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Here  P  denotes  the  principal  value  of  the  integral.  A  detailed 
derivation  of  equation  (4)  can  be  found  in  Cole  (1980).  This  result 
requires  that  G  has  a  specific  asymptotic  behavior  at  its  singularity. 
The  function  G  used  in  the  Kirchhoff  formulation  here  meets  this 
requirement;  Specifically 


GfH.Io.t-V 
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Substitution  of  this  function  into  equation  (2)  gives  a  familiar  optics 
formula  (Born  and  Wolft  1964) 
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where  .r  “  |r-r0l,  the  distance  from  receiver  to  surface,  -  Vu.  jn  and 
■fe  "  V|t-I0|.JL*  The  square  brackets  denote  the  values  of  the  functions 
on  3V  at  the  time  t  -  li~l0l/c. 

The  integrals  (2),  (4),  and  (6)  are  formally  exact  and  are  a 
mathematical  representation  of  Huygen's  principle;  that  is,  a 
disturbance  at  a  receiver  point  is  a  superposition  of  secondary  waves 
proceeding  from  a  surface  existing  between  that  point  and  the  source. 
Diffraction  phenomena  arise  from  the  mutual  Interference  of  these 
secondary  disturbances.  However,  one  needs  the  value  of  the  potential 
and  its  normal  derivative  on  the  surface  to  calculate  u(r^,t).  Equation 
(4)  may  be  solved  for  u(£,t)  or  Vp*n  on  the  surface  subject  to  some 
constraints  imposed  by  boundary  conditions  (e.g.  continuity  of  u(£,t) 
or  V  u»n  across  the  boundary).  This  approach  is  taken  by  Mitaner  (1967) 
who  sets  V  u»  n  equal  to  aero  and  solves  for  u(t,t).  However  this 
approach  may  be  costly  for  high  frequency  scattering. 

Alternatively  one  may  estimate  the  values  on  the  surface  by 
Invoking  an  approximation.  This  approach  is  used  in  this  paper. 
Assuming  a  point  source,  the  boundary  values  on  the  surface  are  taken  to 
be 

u(r,t)  .  (f  (z-  r)/roj(l  +H) 

IS- ■£(>-  ■)(-*(-  ^o2-  £ 

Here  r0  is  the  distance  from  the  source  to  the  surface,  R  is  the 
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approximate  plane-wave  flat  interface  reflection  coefficient,  and  f(t) 
and  f'(t)  are  the  source  time  function  and  its  derivative,  respectively. 
The  reflection  coefficient  will  depend  on  incidence  angle.  This 
approximation  Is  variously  called  the  Kirchhoff,  physical  optics,  or  the 
tangent  plane  hypothesis  and  is  widely  used  by  workers  in 
electromagnetic  scattering  investigations  (Davies, 1954).  It  assumes 
that  the  incident  pulse  is  sufficiently  high  frequency  so  that  locally 
the  amplitude  decay  is  described  by  both  geometric  ray  theory  and 
plane-wave  reflection  coefficients.  Therefore  every  point  on  the 
surface  reflects  the  incident  pulse  as  though  there  were  an  Infinite 
plane  tangent  to  the  surface  at  that  point.  The  values  of  the  potential 
and  its  normal  derivative  at  a  point  is  Independent  of  the  boundary 
values  at  other  points.  Hence  the  effects  of  multiple  scattering  and 
diffractions  along  the  surface  are  neglected. 

Upon  substituting  the  values  (7)  and  (8)  for  u  and  ^  one  obtains 
for  the  reflected  potential 


This  equatic  >  is  similar  to  those  derived  by  A.V.Trorey  (1970, 
1977),  F.J.  Hilterman  (1970,  1975),  and  Berryhill  (1977).  These 
authors  have  derived  convolutional  forms  of  the  Kirchhoff  integral  which 
make  computation  rapid.  Hilterman  has  verified  his  results  with  small 
scale  experimental  modeling  of  a  point  source  in  air  impinging  on  rigid 
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anticlines,  synclines,  and  normal  faults.  The  agreement  between  his 
numerical  calculations  and  experiments  is,  in  general,  excellent. 
However,  these  analytical  forms  of  the  solution  place  severe 
restrictions  on  either  the  source-receiver  geometry  or  the  surface 
geometry. 

The  method  presented  in  this  paper  differs  from  Trorey  and 
Hllterman  in  that  the  source  and  receiver  are  allowed  to  be  at  separate 
locations,  the  surface  geometry  is  arbitrary,  and  the  integral  is 
approximated  by  the  following  expression 


where 

,<».  i—  Us.  ♦  i_  31 

k  rr  2  r  r2  s» 

o  o 
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An  Important  part  of  the  procedure  is  the  discretization  of  the  surface. 
The  rough  surface  is  specified  by  a  function  z(x^,y^)  where  (x^y^)  is  a 
location  on  a  horizontal  grid  of  regularly  spaced  points  separated  by  a 
distance  Ax  and  Ay  (see  Fig.  1).  From  this  information,  one  can 
readily  calculate  and  Then  the  following  formulae  are  used  to 

calculate  AS^  and  : 
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Here  x8,  yg,  and  zg  are  the  station  coordinates  and  xQ,  y0,  and  zQ  are 
the  source  coordinates.  These  values  are  calculated  for  a  finite 
surface.  However  the  Klrchhoff  integral  is  formally  stated  for  a  closed 
surface.  This  problem  is  circumvented  by  integrating  over  a  closed 
surface  consisting  of  the  part  of  the  surface  one  is  Interested  in  (S  ) 
and  a  portion  of  a  sphere  of  large  radius  (Sg)  (see  Fig.  2).  One  then 
argues  that  the  contributions  from  the  sphere  arrive  at  the  receiver  at 
times  later  than  those  of  Interest  and  neglects  them. 

For  the  time  function  f(t)  we  have  chosen  a  ramp  function.  This 
choice  circumvents  the  problems  of  numerically  simulating  a  delta 
function.  So,  on  each  element  we  add  the  sum  of  a  ramp  function 
multiplied  by  and  a  step  function  multiplied  by  appropriately 
in  time.  That  is,  each  element  is  Illuminated  and  contributes  to  the 
total  response  at  a  time  r  »  (r+r0)/c.  This  two-way  travel  time  is 
calculated  and  the  responses  from  all  the  elements  are  summed 
cumulatively  in  order  of  increasing  t  as  displayed  in  Fig.  3. 

For  problems  presented  in  this  paper,  the  numerical  ramp  response 
is  convolved  with  the  analytical  third  derivative  of  a  modified  Haskell 
explosion  source  function  ,  specifically 


(14) 


*(t)  -  Fo  [l  -  e"kt(l  ♦  kt  +  (kt)2/2  -  B(kt)3)] 
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Here  1'0  is  the  source  strength,  k  scales  inversely  as  the  source 
strength  and  B  is  the  overshoot  constant.  This  convolution  is 
mathematically  equivalent  to  the  first  derivative  of  the  reflected 
potential  caused  by  a  source  described  by  (14).  From  trial  and  error, 
we  have  determined  that  such  a  convolution  eliminates  the  splkiness 
introduced  by  simple  differencing  and,  therefore,  is  the  preferred 
method  for  differentiation  in  the  calculations  presented  in  this  paper. 

We  have  carried  out  experiments  to  determine  the  grid  size  required 
to  produce  a  smooth  seismogram.  As  an  example  we  show,  in  figure  4,  the 
variation  of  the  waveform  and  maximum  amplitude  of  a  seismogram  as  a 
function  of  grid  size  for  a  sample  reflection  problem.  The  reflecting 
flat  surface  is  specified  for  seismograms  A,  B,  and  C  by  grid  areas  of  a 
wavelength  of  4  km.  Hence  for  seismograms  A,  B,  and  C  the  number  of 
grids  per  wavelength  is  11.4,  8  and  4  respectively.  Seismogram  C  shows 
that  the  coarse  discrimination  of  the  surface  has  introduced  high 
frequency  noise  which  degenerates  the  waveform  and  the  maximum 
amplitude.  For  each  problem  we  calculate,  the  grid  size  is  selected  by 
trial  and  error.  The  grid  is  made  progressively  finer  until  the 
solution  is  unvarying  as  Bhown  in  figure  4. 

Applications 

The  development  of  any  numerical  procedure  necessitates  a 
comparison  with  known  exact  solutions.  To  see  that  the  integral  is 
being  calculated  correctly,  the  results  of  this  code  are  compared  with 
an  analytical  formula  developed  by  Hllterman  (1975).  For  the  source  and 
the  receiver  together  above  an  arbitrary  rigid  surface,  Hllterman 
reduces  the  Kirchhoff  integral  to  the  following  form: 
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(15) 


Here  *  denotes  convolution.  ^  is  the  Increase  of  solid  angle  with  a 
vertex  at  the  source-receiver  point  subtended  by  the  surface  as  a 
function  of  time  (Hilterman,  1975).  We  choose  to  calculate  the  response 
to  a  point  source  in  a  fluid  overlying  a  rigid  hemisphere  imbedded  in  a 
infinite  rigid  plane.  (see  Fig.  5  for  the  geometry)  For  such  a 
geometry  we  have  determined  ^  to  be 
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where  tq  ■  2( |&|  -a)/c,  the  minimum  two-way  travel  time,  and  Tj  ■  2( |£| ^ 
+  a2)l/2/c.  H(t)  is  the  assumed  source  time  function  and  is  the  unit 
step  function.  Using  these  results,  we  calculate  the  following  two 
terms  $A  and  analytically  and  numerically  to  check  on  the  accuracy  of 
the  Integral  approximation. 
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The  comparison  is  shown  in  Figure  5.  In  this  calculation,  the 
source-receiver  point  is  20  km  above  the  center  of  the  hemisphere  with  a 
radius  of  5  km.  The  velocity  of  the  medium  is  6  km/sec.  The  agreement 
is  good  for  tq  +  IS  seconds.  The  results  differ  because  the  integral  is 
calculated  numerically  over  a  finite  surface.  The  conclusion  from  this 
experiment  is  that  the  numerical  evaluation  of  the  integral  is  adequate. 

We  further  test  the  code  by  comparing  the  Kirchhoff  solutions  with 
analytical  high  frequency  solutions  which  satisfy  the  given  boundary 
conditions.  Again  we  choose  to  calculate  the  potentials  caused  by  a 
point  source  impinging  on  a  rigid  sphere;  however  we  allow  the  source 
and  receiver  to  separate.  The  Kirchhoff  results  are  compared  -fo  those 
from  a  first  motion  approximation  of  an  aymptotic  form  of  a  solution 
obtained  by  Gilbert  and  Helmberger  (1972).  They  solve  the  problem  of 
the  reflection  of  an  SH  pulse  from  a  fixed  and  rigid  sphere.  Figure  6 
illustrates  the  geometry  and  parameters  used  in  this  problem.  The 
displacement  as  a  function  of  spherical  polar  coordinates  (r,0)  obeys 
the  following  equation. 

ot  2*r2sin6 


Here  p  is  the  rigidity  of  the  medium,  p  is  the  density,  and  c  is  the 
wavespeed  (c  -(/P>.  The  Gilbert  and  Helmberger  asymptotic  solution  for 
the  reflected  pulse  is 
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f(s)  are  the  Laplace  transform  of  u(r,8,t)  and  f(t).  The  variables 
Y  and  ♦  and  the  path  of  Integration  are  defined  in  Gilbert  and 
Helmberger  (1972).  After  performing  the  first  motion  approximation  we 
obtain  for  the  reflected  pulse 
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and  for  the  direct  pulse 
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where  R'  is  the  distance  between  the  source  and  the  receiver. 

If  the  source  and  receiver  are  together,  it  can  be  shown  that 
solution  (22)  and  a  far-field  first  motion  approximation  to  the 
Kirchhoff  solution  are  equivalent.  First  let  us  define  a  geometric 
spreading  factor  S  such  that  equation  (22)  may  be  rewritten  as 
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Now  we  examine  the  Kirchhoff  solution  when  the  source  and  receiver  are 
together.  Following  an  approach  developed  by  Hilterman  the  first  term 
in  equation  (15)  Is  discarded  as  a  near-field  term. 
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Then  a  first  motion  approximation  is  made. 
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Then 
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From  equation  (16)  we  find  that 
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If  we  now  take  the  limiting  value  of  S  as  R  R0  and  6  ♦  0  In  equation 
(22),  the  result  is 

8  “  2(R  +  a)R 
o  o 


(29) 


Hence,  the  two  far-field  high  frequency  solutions  are  formally 
equivalent  when  the  source  and  receiver  are  together.  A  similar  result 
is  obtained  by  Hllterman  (1975)  for  a  rigid  planar  surface. 

The  Kirchhoff  solutions  are  also  calculated  when  the  source  and 
receiver  are  separated  and  the  maximum  amplitudes  of  the  synthetics  are 
compared  with  those  amplitudes  predicted  by  equation  (22).  The  time 
history  of  the  input  source  for  this  problem  is  described  by  the  third 
derivative  of  equation  (14)  with  the  overshoot  constant  B  ■  2  and  k  * 
10.  The  medium  has  a  shear  wave  velocity  of  5  km/sec;  thus  the 
wavelength  of  the  input  source  is  approximately  4  km.  The  grldlength 
used  in  these  calculations  is  .1  km,  making  the  number  of  grids  per 
wavelength  equal  to  40.  The  total  grid  area  needed  to  describe  the 
surface  of  the  sphere  is  400  km  .  The  grldlength  was  selected  to  give 


an  extremely  fine  sampling  of  the  surface  so  that  we  may  investigate  the 
effects  of  a  wide  range  of  pulse  widths  as  Input  time  histories.  Six 
ramp  responses  for  this  problem  required  595. 8  seconds  of  CPU  time  on  a 
PRIME750.  Table  1  shows  the  parameters  used  in  tt  numerical 

experiments  and  the  numerical  and  theoretical  amplitudes.  The  results 
compare  favorably;  the  accuracy  of  the  Kirchhoff  solutions  is  better 
than  or  equal  to  1  Z. 

The  above  two  experiments  indicate  that  the  Kirchhoff  code 
correctly  predicts  reflections  from  curved  surfaces  with  large 
reflection  coefficients  and  far-field  receivers.  Similiar  efforts  have 
been  carried  out  by  workers  in  the  field  of  electromagnetic  scattering. 
Jlracek  (1972)  computes  the  amplitudes  of  electromagnetic  waves  caused 
by  an  incident  transverse  electric  plane  wave  impinging  on  a  perfectly 
conducting  two-dimensional  sinusoidal  surface.  He  compares  results 
obtained  from  a  Rayleigh-FFT  method,  an  integral  equation  solver,  and 
the  Kirchhoff  method.  The  most  obvious  failure  of  the  Kirchhoff 
technique  to  predict  correct  amplitudes  occurs  when  the  incident  angle 
is  past  critical  angle.  This  result  is  not  surprising  in  light  of 
assumptions  made  in  estimating  the  boundary  values  on  the  surface. 
However,  for  angles  less  than  critical,  the  Kirchhoff  code  is  adequate 
and  inexpensive  for  problems  Involving  three-dimensional  rough  surfaces. 

Ve  can  gain  further  insight  into  the  usefulness  of  this  formalism 
by  comparison  of  the  Kirchhoff  solutions  with  optical  solutions  to 
problems  of  geophysical  Interest.  First,  the  technique  is  applied  to 
the  calculation  of  reflections  from  a  mountain  with  a  buried  source.  In 
the  second  application,  reflections  from  a  plane  where  the  reflection 
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coefficient  varies  as  a  function  of  position  on  the  surface  are 
computed.  In  both  calculations,  particular  attention  will  be  paid  to 
those  propagation  paths  where  classical  ray  theory  fails. 

The  first  application  of  the  code  is  the  calculation  of  the 
reflected  potentials  from  an  isotropic  source  underneath  an  Idealized 
mountain  (see  Fig.  7).  The  topography  csf  the  mountain  is  calculated 
from  this  formula  where  1  is  the  height  of  the  surface. 

Here  C,  the  maximum  height,  is  5  km,  and  W,  the  width,  is  33.33  km.  The 
acoustic  reflection  coefficient  is  -1  everywhere  on  the  surface.  The 
topography  is  specified  on  150  x  150  km  grid  with  each  element  of  the 
grid  being  .5  km  long. 

Since  the  angle  between  the  normal  of  the  surface  and  the  Incident 
source  ray  is  calculated  by  the  code,  it  is  simple  to  plot  the  path  of 
the  reflected  rays.  These  rays  are  traced  for  two  depths  below  the 
baseline  of  the  free  surface.  In  the  first  plot.  Figure  7,  we  can  see 
the  rays  from  a  source  at  10  km  which  reflect  off  the  free  surface  and 
travel  to  a  depth  of  50  km.  This  figure  shows  the  position  of  the  ray 
caustic,  focil,  and  the  shadow  zone  caused  by  the  convex  shape  of  the 
mountain.  These  features  will  Influence  the  waveforms  considerably. 

In  Figure  8  the  rays  are  traced  to  a  depth  of  1000  km.  The 
Klrchhoff  responses  are  calculated  at  this  depth  at  the  marked  positions 
which  vary  from  0  kilometers  to  750  kilometers  horizontally.  Upon 
closer  inspection  of  Figure  9,  one  can  see  slight  asymmetries  in  the 
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location  of  the  rays  with  respect  to  the  position  at  0  kilometers. 
These  asymmetries  are  caused  by  the  discretization  of  surface  of  the 
mountain.  The  error  in  the  value  of  the  competed  normal  derivatives 
introduces  10  km  of  uncertainty  into  the  location  of  the  rays  at  this 
depth. 

The  calculated  reflected  responses  are  shown  in  Figure  9.  These 
pulses  are  convolutions  of  the  ramp  response  with  the  Haskell  function 
with  the  parameter  B-0  and  k«25.  Hence  the  number  of  grids  per 
wavelength  is  10.  As  the  horizontal  distance  of  the  receiver  changes, 
we  see  systematic  waveform  variations  which  can  be  interpreted  In  terms 
of  rays  Interacting  with  caustics.  In  the  ranges  of  0,  50,  100,  and  150 
kilometers  the  synthetics  have  complicated  pulse  shapes  caused  by  the 
Interference  of  three  families  of  rays.  The  first  arrival  is  a  simple 
pulse  with  a  *  phase-shift  which  is  a  consequence  of  the  reflection  off 
the  free  surface.  The  second  pulse  is  a  reflected  ray  with  a  path  which 
is  tangent  to  the  caustic  formed  by  the  mountain.  This  path  results  in 
a  i  +  ir/2  phase-shift  of  the  pulse.  The  third  arrival  reflects  off  the 
mountain  and  travels  through  the  geometric  focus  caused  by  the  mountain; 
thus  the  phase  shift  of  this  arrival  1b  it  +  it  ,  The  maximum  amplitude 
of  these  four  distances  is  controlled  by  the  interference  of  these  rays. 
Clearly  the  high  amplitude  and  the  simple  pulse  of  the  first  synthetic 
at  0  kilometers  is  a  result  of  the  constructive  Interference  of  the 
first  two  rays.  Fast  150  kilometers,  the  latter  two  arrivals  arrive 
closely  in  time  and  their  Interference  controls  the  amplitude  and 
frequency  content  of  the  second  pulse  on  the  record.  From  Figure  8  ,  it 
is  clear  that  a  ray  interpretation  of  pulses  on  records  past  400 
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kilometers  Is  no  longer  valid.  Ray  theory  predicts  only  one  reflected 
pulse  because  the  mountain  creates  a  shadow  zone;  yet  one  sees  two 
distinct  pulses  predicted  by  the  Klrchhoff  method.  The  second 
phase-shifted  pulse  decreases  in  amplitude  and  frequency  content.  As 
the  horizontal  distance  of  the  receivers  increases,  the  amplitude  of  the 
first  reflection  becomes  the  maximum  amplitude  of  the  record.  If  one 
calculates  the  maximum  amplitudes  of  reflections  off  a  plane  for  the 
same  source-  receiver  geometry,  one  can  see  that  the  two  sets  of 
amplitudes  merge.  This  behavior  is  shown  in  Figure  10  where  the 
amplitudes  as  a  function  of  horizontal  distance  for  the  two  geometries 
have  been  calculated.  The  solid  line  shows  the  decay  of  amplitudes 
calculated  for  a  planar  surface.  The  triangles  are  amplitudes 
calculated  for  a  mountain  with  a  height  of  2  kilometers  and  a  width  of 
10  kilometers.  The  two  sets  of  values  coincide  past  800  kilometers. 

The  Klrchhoff  results  for  this  experiment  are  gratifying  because 
one  does  not  expect  infinite  amplitudes  or  abrupt  shadow  zones  predicted 
by  optics  in  real  physical  systems.  This  experiment  also  demonstrates 
that  this  technique  produces  the  requisite  phase  shifts  in  an  extremely 
simple  manner  unlike  existing  ray  tracing  techniques  which  must  track 
the  behavior  of  a  ray  tube  along  the  propagation  path. 

The  second  application  of  the  code  is  the  calculation  of 
reflections  off  an  acoustic  planar  free  surface  where  the  reflection 
coefficients  are  allowed  to  vary  as  a  function  of  position  on  the 
surface.  These  calculations  demonstrate  the  flexibility  of  the  code  and 
again  emphasize  the  differences  between  the  Klrchhoff  solution  and 
optics.  (The  wavespeed  of  the  medium  is  6km/sec  for  all  the  following 
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calculations). 

Initially  one  assumes  that  the  reflection  coefficient  is  zero  for 
elements  of  the  plane  within  a  circular  aperture  of  radius  R  and  is  -1 
for  elements  outside  this  aperture.  The  source  la  directly  underneath 
the  center  of  the  hole.  From  ray  theory  one  expects  that  no  reflected 
energy  will  arrive  at  a  receiver  directly  underneath  the  source.  Yet 
one  calculates  non-zero  amplitudes  for  both  long  and  short  period  WWSSN 
seismograms  from  the  Klrchhoff  code.  These  seismograms  are  displayed  in 
Figure  11  as  a  function  of  the  radius  of  the  aperture  for  a  receiver 
1000  kilometers  below  the  surface.  Only  the  reflected  pP  phase  has  been 
calculated  and  convolved  with  WWSSN  instruments. 

This  pulse  is  systematically  delayed  as  the  radius  of  the  hole 
increases  from  1  to  5  kilometers.  There  is  no  change  in  the  waveforms. 
Only  the  amplitudes  of  both  sets  of  seismograms  decrease.  However  the 
amplitude  of  the  seismograms  for  an  aperture  with  a  five  kilometer 
radius  is  more  than  half  the  magnitude  of  the  amplitude  of  a  pP  phase 
reflected  from  a  free  surface  without  a  hole.  Clearly,  then,  ray  theory 
is  not  a  good  approximation  to  the  solution  to  this  problem. 

In  addition,  ray  theory  falls  to  predict  any  dependence  of  the 
reflected  amplitudes  on  frequency.  Intuitively  one  expects,  for  an 
aperture  problem,  that  the  higher  frequencies  of  a  broad-band  signal 
will  be  reduced  relative  to  the  lower  frequencies  after  reflection. 
This  hypothesis  is  tested  by  calculating  the  reflected  responses  from 
sources  of  differing  frequency  content.  In  the  following  calculations 
the  parameter  B  of  the  modified  Haskell  source  representation  equals 
zero;  however  k  varies  from  5  to  25.  An  Increase  in  k  broadens  the 
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bandwidth  of  the  incident  signal  (von  Seggern  and  Blandford,  1972).  One 
computes  two  responses  for  a  given  source  pulse.  The  first  response  is 
a  reflection  off  the  plane  with  a  hole  and  the  second  i6  a  reflection 
off  the  plane  without  a  hole.  The  amplitude  of  the  latter  response  has 
no  frequency  dependence;  hence  if  the  reflection  from  a  hole  has  no 
frequency  dependence,  one  predicts  that  the  ratio  of  the  amplitudes  of 
the  two  reflections  will  be  constant  as  a  function  of  the  parameter  K. 
However,  if  there  is  a  frequency  dependence,  the  ratios  should  vary 
systematically. 

From  numerical  experiments  one  confirms  the  frequency  dependence  of 
the  reflected  amplitudes.  This  result  is  shown  in  Figure  12. 
Specifically,  the  amplitudes  of  the  reflections  from  the  aperture  are 
always  smaller  than  the  amplitudes  of  planar  reflections.  Also  the 
ratio  of  the  two  responses  decreases  when  B  decreases  if  the  receiver  is 
located  at  position  2.  This  behavior  is  displayed  for  apertures  with 
three  radii,  2,  3,  and  A  tat. 

However,  this  behavior  does  not  occur  If  the  receiver  is  located  at 
position  1,  1000  tat  directly  below  the  source.  The  ratios  are 
approximately  constant  as  B  decreases.  This  observation  suggests  that, 
here,  the  reflected  amplitudes  from  the  aperture  do  not  depend  on  the 
bandwidth  of  the  signal.  Although  this  result  is  not  intuitive,  it  is 
typical  of  analytical  solutions  to  Fraunhofer  diffraction  from  apertures 
in  an  opaque  screen  (Born  and  Volf,  1964).  For  example,  the  solution  of 
the  intensity  of  light  transmitted  through  a  rectangular  aperture  has 
the  functional  form  of 

sin  (Atux)  9  sin  (Buy) 

Aux  Buy 


(31) 
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where  A  and  B  are  geometric  constants  and  x  and  y  are  the  rectangular 
coordinates  of  the  position  of  the  receiver.  The  limit  of  the  above 
function  as  x  and  y  approach  zero  is  1  and  is  independent  of  the  value 
of  o. 

Such  experiments,  which  vary  the  reflection  coefficients  on  the 
free  surface,  may  be  applicable  to  the  analysis  of  the  effects  of 
spallation  generated  by  nuclear  blasts  on  teleseismic  P-wave 
reflections.  Spall  is  the  physical  separation  of  near  surface  layers 
during  the  explosion.  Material  above  the  bomb  is  either  ejected  or 
returns  to  produce  an  Impact  signal  on  near-field  Instruments.  This 
non-linear  and  non-elastic  behavior  of  the  material  surrounding  the 
source  may  result  in  amplitude  and  travel  time  anomalies  of  reflected  sP 
and  pP  phases. 

The  model  used  to  simulate  spall  is  one  where  the  reflection 
coefficient  is  a  cosine  taper;  that  is 

R  -  (c°»(Jf)-l)  *<K 


R  •  -A  X  >  R  (33) 

Here  x  is  distance  from  the  source  epicenter  on  the  free  surface.  One 
chooses  this  behavior  of  the  reflection  coefficient  to  simulate  material 
reflecting  more  energy  as  the  distance  from  the  source  increases.  The 
model  Introduces  complications  Into  the  short  period  waveforms  but  only 
broadens  the  long  period  waveforms.  This  effect  and  the  source-receiver 
geometry  is  illustrated  in  Figure  13.  The  geometry  is  the  same  as  used 
in  the  aperture  calculations.  Unlike  the  first  model,  this  model  causes 
the  amplitudes  of  both  the  long  and  short  period  reflections  to  decay 


quite  rapidly.  The  amplitude  decay  is  greater  for  the  short  period 
reflections  than  for  the  long  period  reflections.  Hence  the  long  period 
energy  is  insensitive  to  the  perturbation  of  the  reflection  coefficients 
relative  to  the  short  period  energy. 

The  source-receiver  geometry  is  changed  for  this  particular  model 
to  test  the  hypothesis  that  asymmetries  of  spalling  with  respect  to  the 
source  location  can  Introduce  observable  azimuthal  variations  of 
amplitudes  and  waveforms  of  teleseismic  records  of  nuclear  blasts.  Such 
variations  have  been  documented  for  teleseismic  recordings  of  Nevada 
Test  Site  blasts  (Helmberger  and  Hadley,  1981).  In  addition, 
photographs  of  collapse  craters  from  NTS  blasts  suggest  that  processes 
like  spalling  and  subsidence  occur  along  pre-existing  planes  of  weakness 
which  are  not  symmetrical  with  respect  to  the  emplacement  hole  (Springer 
and  Kinnaman,  1971).  Figure  14  shows  the  results  for  stations  at  three 
azimuths.  The  source  is  placed  2  kilometers  to  the  right  of  the  center 
of  the  spall  aperture  and  1  kilometer  below  the  free  surface.  The 
receivers  are  all  at  horizontal  distance  which  corresponds  to  a  take-off 
angle  of  20°  for  the  direct  P  wave.  One  sees  azimuthal  variations  of 
waveform  and  amplitudes  for  both  long  and  short  period  reflections.  The 
amplitude  variations  are  not  large  but  the  waveform  changes  are  dramatic 
for  short  period  records. 

This  model  is  crude  and  consequently  does  not  demonstrate  that 
spalling  affects  teleseismic  reflections.  However,  Shumway  and 
Blandford  (1980'*  report  observing  a  systematic  delay  in  arrival  times  of 
pP  phases  from  explosions.  The  simple  aperture  experiment  provides  an 
explanation  for  that  delay.  In  addition  the  Kirchhoff  technique  allows 
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one  to  specifiy  more  realistic  dynamical  information  on  the  free  surface 
and  calculate  more  realistic  models  in  a  straightforward  manner. 

Conclusions 

A  numerical  procedure  has  been  presented  for  the  evaluation  of  the 
Kirchhoff-  Helmholtz  integral  assuming  the  tangent  plane  hypothesis. 
The  method  is  a  high  frequency  one  and  produces  results  which  compare 
well  with  existing  asymptotic  first  motion  solutions.  The  technique  has 
been  applied  to  two  problems  and  compared  to  classical  ray  theory 
results.  First,  the  reflections  off  an  idealized  mountain  are 
calculated  and  have  phase-shifts  consistent  with  those  predicted  by 
optics;  however,  the  amplitudes  at  triplications  are  finite  unlike  the 
classical  ray  result.  In  addition  diffracted  pulses  are  produced  in  the 
shadow  zones.  The  second  application  is  the  calculation  of  reflections 
where  the  reflection  coefficients  vary  as  a  function  of  position.  For  a 
hole  in  the  free  surface,  the  Kirchhoff  method  produces  reflections 
where  ray  theory  predicts  no  reflections.  The  method  also  produces 
amplitudes  which  are  frequency  dependent.  The  results  are  applied  and 
extended  to  model  the  effects  of  spallation  on  teleseismic  reflections. 
Travel  time  delays  and  amplitude  anomalies  are  predicted.  These 
anomalies  are  consistent  with  observations  although  the  observations  are 
not  modeled. 

In  conclusion,  the  method  has  a  broad  range  of  applications.  The 
method  is  inexpensive  to  run  for  modeling  two  and  three  dimensional 
rough  surfaces.  Although  the  method  is  appropriate  for  narrow  angles  of 
reflections  and  acoustic  reflections,  its  range  of  applicability  can  be 
extended  by  introducing  new  time  functions  on  the  boundary.  The  code 
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can  also  be  coupled  with  existing  propagational  techniques  such  a6 
ray-tracing,  Cagniard-de  Hoop  methods,  or  full  wave  theory.  This 
coupling  will  enable  one  to  handle  more  complicated  and  relevant 
selsmological  problems. 
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FIGURE  CAPTIONS 

The  discretization  of  the  surface  over  which  the  Klrchhoff 
Integral  is  calculated. 

The  closed  surface  of  Integration  which  is  composed  of  Sa  the 
surface  over  which  one  carries  out  the  numerical 
integration,  and  Sg  ,  the  remainder  which  does  not 
contribute  to  the  signal  in  the  time  window  of  interest. 

The  summation  of  the  response  from  each  element  to  obtain  the 
total  response. 

This  figure  shows  the  change  of  response  as  the  element  size 
in  creases.  The  model  used  here  is  a  free  surface  with  a 
hole.  The  source  and  receiver  are  directly  below  this  hole. 
The  response  is  convolved  with  a  modified  Haskell  source  and 
a  short  period  WWSSN  Instrument  to  produce  the  above 
seismograms. 

The  geometry  of  the  point  source  problem  is  shown  at  the  top 
of  figure.  Here  a  is  the  radius  of  the  hemisphere.  r_  is 
the  vector  from  the  sour ce -receiver  point  to  an  arbitrary 
position  on  the  surface.  R  is  the  vector  extending  from  the 
center  of  the  hemisphere  to  the  source-receiver  point. 
Below  are  the  two  parts  of  the  solution  and  .  The 
dotted  lines  are  the  values  computed  by  the  numerical 
Integration.  The  solid  lines  are  the  values  computed  by  the 
analytical  Hilterman  solution. 
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6.  Geometry  of  the  spherical  problem. 

7.  The  rays  which  reflect  off  a  mountain  described  by  the 

equation  (30).  The  source  is  10  km  below  the  baseline 
Indicated  by  the  dashed  line.  Also  shown  are  the  caustics 
formed  by  such  a  mountain  and  the  geometric  focus. 

8.  The  rays  which  reflect  off  the  mountain  are  traced  to  a  depth 

of  1000  km. 

9.  The  responses  convolved  with  a  modified  Haskell  source  for 

receivers  all  at  a  depth  of  1000  km  and  at  horizontal 
distances  which  vary  from  0  to  750  km  away  from  the  center 
position. 

10.  The  maximum  amplitude  of  the  waveforms  calculated  for  two 

free  surfaces  is  shown  as  a  function  of  horizontal  distance 
away  from  the  center.  The  receivers  are  at  a  depth  of  1000 
km.  The  solid  line  is  the  value  of  the  maximum  amplitude  of 
reflections  off  a  plane  from  a  source  10  km  below  the  free 
surface.  The  triangles  are  values  for  a  free  surface  with  a 
mountain  of  height  2  km  and  width  10  km  for  the  same 
source-receiver  geometry. 

11.  Long  (dashed  lines)  and  short  (solid  lines)  period  VWSSN 

seismograms  calculated  for  the  geometry  shown  at  the  top  of 
the  figure.  Here  the  radius  of  the  hole  varies  from  1  to  5 
km. 

12.  Ratio  of  maximum  amplitude  of  reflection  from  a  surface  with 

and  without  a  hole.  The  open  circles  are  ratios  measured  at 
a  receiver  directly  underneath  the  source  and  the  center  of 
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the  aperture.  The  triangles  are  ratios  measured  at  a 
receiver  68  km  away  from  the  center  position.  Both 
receivers  are  1000  km  be  low  the  free  surface.  The  radius 
of  the  hole  R  is  2,3,  and  4  km.  The  ratios  are  plotted  for 
these  radii  as  a  function  of  k.  (a  measure  of  the 
bandwidth) 

13.  Long  and  short  period  seismograms  for  a  cosine  tapered 

reflection  coefficient. 

14.  Long  and  short  period  seismograms  from  receivers  located  at 

positions  A  ,  B,  and  C. 
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Parameters  used  in  calculation  of 
Response  from  a  Sphere  with  radius 
a=  10km 
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WWSSN  SYNTHETICS 
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